NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


THESIS 

NUMERICAL  OPTIMIZATION  ALGORITHM 

FOR  ENGINEERING  PROBLEMS 

USING  MICROCOMPUTER 

by 

Dong  Soo,  Kim 

September  1984 

Thesis  Advisor:                  G.  N.  Vand( 

2rplaats 

Approved  for  public  release;  distribution  unlimited. 


1221 J k^ 


Unclassified 


SECURITY  CLASSIFICATION   OF   THIS  PAGE  C»h*n  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


1.  REPORT  NUMBER 


2.  GOVT  ACCESSION  NO 


3.  RECIPIENT'S  CATALOG  NUMBER 


4.     TITLE  ("and  Subt/r/e; 


Numerical  Optimization  Algorithm  for 
Engineering  Problems  Using  Microcomputer 


5.     TYPE   OF    REPORT   4    PERIOD   COVERED 

Master's  Thesis; 
September  1984 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHORl'a; 


8.  CONTRACT  OR  GRANT  NUMBERfs; 


>Dong    Soo ,    Kim 


9.     PERFORMING  ORGANIZATION   NAME   AND   ADDRESS 

Naval  Postgraduate  School 
Monterey,  California   93943 


10.     PROGRAM   ELEMENT,  PROJECT,   TASK 
AREA  ft   WORK  UNIT   NUMBERS 


tl.     CONTROLLING   OFFICE   NAME   AND   ADDRESS 

Naval  Postgraduate  School 
Monterey,  California   93943 


12.      REPORT   DATE 

September   1984 


13.     NUMBER  OF   PAGES 

57 


U.     MONITORING  AGENCY  NAME  4    ADDR ESSf// dy/Zoronf  Itom  Controlling  Oltlca) 


15.     SECURITY   CLASS,   (ol  thia  report) 


Unclassified 


15«.     DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


16.     DISTRIBUTION   STATEMENT  (of  this  Report) 


Approved  for  public  release;  distribution  unlimited, 


17.     DISTRIBUTION   STATEMENT  (ot  the  abstract  entered  In  Block  20,  II  different  from  Report) 


18.     SUPPLEMENTARY  NOTES 


19.     KEY   WORDS  (Continue  on  reverse  aide  If  neceaaary  and  Identity  by  black  number) 

Microcomputer 

Feasible  Direction  Method 


20.      ABSTRACT  (Continue  on  reverse   aide   If  necesaary   and  Identify  by  block  number) 

A  general  purpose  computer  program  is  developed  to  perform 
nonlinear  constrained  optimization  of  engineering  design  problems. 
The  program  is  developed  especially  for  use  on  microcomputers  and 
is  called  Microcomputer  Software  for  Constrained  Optimization 
Problems  (MSCOP) .   It  will  accept  a  nonlinear  objective  function 
and  up  to  50  inequality  constraint  functions  and  up  to  20  bounded 
design  variables. 


DD     1   JAN   73     1473  EDITION  OF    1   NOV  65  IS  OBSOLETE 

S    N  0102-  LF-  0)4-  6601  1 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Sntarad) 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Wh»n  Data  Bnffd) 


MSCOP  employs  the  method  of  feasible  directions.   Although 
developed  for  microcomputers,  for  speed  of  development,  the 
MSCOP  was  implemented  on  an  IBM  3033  using  standard  basic 
language,  Waterloo  BASIC  Version  2.0.   It  is  directly  trans- 
portable to  a  variety  of  microcomputers. 

Typical  applications  of  MSCOP  program  are  in  the  design  of 
machine  components  and  simple  beam  and  truss  structures.   Solutions 
to  three  sample  problem.s  are  given. 


S    N   0102-  LF-  014-  6601 

2  Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE(T»Ti«n  Dalm  Enfrmd) 


Approved  for  public  release;  distribution  unlimited, 


Numerical  Optimization  Algorithm 

for  Engineering  Problems 

Using  Micro-computer 


by 


Dong  Soo,  Kim 
Major,  Republic  of  Korea  Army 
B.S.,  Korea  Military  Academy,  1976 
B.E.,  Seoul  National  University,  1980 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  1984 


DUOLt^_^-;-^^£  SCHOOL 


ABSTEACT 

A  general  purpose  computer  program  is  developed  to 
perform  nonlinear  constrained  Oiptimization  of  engineering 
design  problems.  The  program  is  developed  especially  for 
use  on  iricrocoraputers  and  is  called  Microcomputer  Software 
for  Constrained  Optimization  Problems  (msco?) .  It  will 
accept  a  nonlinear  objective  function  and  ap  to  50 
inequality  constraint  functions  and  up  to  20  bounded  design 
variables. 

MSCOP  employs  the  method  of  feasible  directions. 
Although  developed  for  microcomputers,  for  speed  of  develop- 
ment, the  MSCOP  was  implemented  on  an  IBM  3033  using  stan- 
dard basic  language,  Waterloo  BASIC  Version  2.0.  It  is 
directly  transportable  to  a  variety  of  microcomputers. 

Typical  applications  of  MSCOP  program  are  in  the  design 
of  machine  components  and  simple  beam  and  truss  structures. 
Solutions  to  three  sairple  problems  are  given. 
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I.  INTIOPnCTION 

A.   PDEPCSE 

This  thesis  describes  the  development  of  a  micr ccomputer 
oriented  program  called  ?1SC0P  (Microcomputer  Software  for 
Constrained  Optimization  Problems)  for  constrained  optimisa- 
tion of  engineering  design  protlems.  Problems  which  can  be 
solved  by  the  MSCCP  are  nonlinear  programming  problems 
arising  in  several  areas  of  machine  and  structural  design, 
such  as  the  minimum  weight  design  of  structures  subject  to 
stress  and  displacement  constraints  [Ref.  1]. 

In  recent  years,  several  powerful  general  purpose  opti- 
mization programs  have  become  available  for  engineering 
design  problems,  e.g.,  COPES/CONMIN  [Ref.  2],  and  ADS-1 
[Ref.  3].  These  programs  can  handle  a  wide  range  of  design 
problems  and  contain  a  variety  of  solution  techniques. 
Also,  several  programs  are  available  that  include  optimiza- 
tion in  an  integrated  analysis  /  design  code,  e.g.,  ACCESS, 
ASO?,  EAI,  PARS,  SAVES,  SPAR,  STARS  and  TSO  [Ref.  U].  All 
of  the  above  optimization  programs  are  written  in  FORTRAN, 
and  are  built  for  use  on  a  mainframe  computer.  Their  use  can 
be  cumbersome,  especially  for  the  occasional  user.  Since 
many  engineers  are  now  using  microcomputers,  there  is  a  need 
to  develop  an  optimization  program  contained  in  a  microcom- 
puter software  package  for  use  on  microcomputers.  This 
thesis  fills  that  need  by  developing  a  compact  program 
written  in  a  standard  BASIC  language  suitable  for  a  wide 
range  of  microcomputers. 
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B.  lEPLEHEHTATTON 

The  nature  of  an  optimization  program  depends  on  the 
computer  and  programming  method  available.  The  MSCO?  soft- 
ware is  designed  for  use  on  a  microcomputer.  However,  for 
the  speed  of  development  and  testing,  MSCOP  was  developed  on 
the  lEM  3033  computer  at  the  F.  P.  Church  Computer  Center  in 
Naval  Postgraduate  School,  and  was  written  in  WEASIC 
(Waterloo  Basic)  Version  2.0. 

To  make  sure  that  the  program  is  easily  portable  to  a 
microcomputer,   only   standard  BASIC  commands   and  functions 

are  used.   For  example,  FOR  I  =  1  TO  ^TDB NEX"^  I,   G0SU3 

etc.,  were  used.  The  commands  and  functions  not  available 
in  all  variations  of  BASIC  are  avoided,  for  example,  TRN(A), 
MAT  (A),  etc. 

nSCOE  provides  design  engineers  with  a  convenient  tool 
for  optinization  of  engineering  design  problems  with  up  to 
20  bounded  design  variables  and  as  many  as  50  inequality 
constraints. 

C.  GENEEAL    OPTIMIZATION    MODEL 

The  general  optinization  problem  to  be  solved  is  of  the 
form  :   Find  the  set  of  design  variables  ^   that  will 

Minimize        F(X)  (1.1) 

Subject  to      G  (X)  <  0  3  =  T/ /HI       (1.2) 

J 

1  u 

X<X<X  i=1, ,    n  (1.3) 

i    -        i   -         i 

where  X  is  referred  to  as  the  vector  of  design  variables. 
F (X)  is  the  objective  function  which  is  to  be  minimized. 
G  (X)  are  inequality  constraint  functions,  and  Xj^  and  X 
are    lower   and      upper    bounds,       respectively,         on      the    design 
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variatles.  Although  these  bounds  or  "side  constraints" 
could  be  included  in  the  inequality  constraint  set  given  by 
Eq(1.2)  ,  it  is  convenient  to  treat  them  separately  because 
of  their  special  structure.  The  objective  function  and 
constraint  functions  may  be  nonlinear,  explicit  or  implicit 
in  X.  However,  they  must  be  continuous  and  should  have 
continuous  first  derivatives. 

In  general  engineering  optimization  problems,  the  objec- 
tive to  be  minimized  is  usually  the  weight  or  volume  of  a 
structure  being  designed  while  the  constraints  gives  limits 
en  compressive  stress,  tensile  stress,  Enler  buckling, 
displacement,  frequencies  (eigenvalues),  etc.  [Eef.  5  : 
p. 264].  Equality  constraints  are  not  included  because  their 
inclusion  complicates  the  solution  techniques  and  because  in 
engineering  situations,  equality  constraints  are  rare. 

Most  optimization  algorithms  require  that  an  initial 
value  of  design  variables  xo  be  specified.  Beginning  from 
these  starting  values,  the  design  is  iteratively  improved. 
The  iterative  procedure  is  given  by 

g+1     q         q 
X     =  X   +   a*  S.  (1.4) 

where  q  is  the  iteration  number,  S  is  a  search  direction 
vector  in  the  design  space,  and  a*  is  a  scalar  parameter 
which  defines  the  amount  of  change  in  X.  At  iteration  q,  it 
is  desirable  to  determine  a  direction  S  which  will  reduce 
the  objective  function  (usable  direction)  without  violating 
the  constraints  (feasible  direction).  After  determining  the 
search  direction,  the  design  variables,  X,  are  updated  by  Eg 
(1.4)  so  that  the  minimum  objective  value  is  found  in  this 
direction.   [Ref.  6]. 

Thus,  it  is  seen  that  nonlinear  optimization  algorithms 
for  the  general  optimization  problem  based  on  Eq{1.4)  can  be 
separated  into  two  parts,  determination  of  search  direction 
and  determination  of  scalar  parameter  a*. 
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D.   OEGAKIZATION  OF  THIS  THESIS 

This  chapter  has  stated  the  purpose  of  the  thesis  and 
has  put  the  general  concept  of  engineering  optimization  into 
a  preliminary  perspective.  Chapter  2  will  describe  the 
essential  aspects  of  the  optimization  algorithm  used  in 
MSCOP  such  as  finding  a  search  direction,  the  one- 
dimensional  search  and  convergence  criteria.  Chapter  3 
describes  program  usage.  In  chapter  4,  there  are  three 
examples  which  are  solved  by  the  MSCOP.  Summary  and  conclu- 
sions are  given  in  chapter  5.  The  program  is  listed  in  the 
appendix. 
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II.  OPTIMIZATION  ALGORITHM 

A.   INTRCDOCTION 

There  are  many  optimization  algorithms  for  constrained 
nonlinear  problems  such  as  generalized  reduced  gradient 
method,  feasible  direction  method,  penalty  function  methods. 
Augmented  Lagrangian  multiplier  method,  and  sequential 
linear  programming.  The  feasible  direction  method  is  chosen 
for  development  in  this  thesis  for  three  main  reasons. 
First  it  progresses  rapidly  to  a  near  optimum  design. 
Second  it  only  requires  gradients  of  objective  and 
constraint  functions  that  are  active  at  any  given  poirt  in 
the  optimization  process  [Ref.  7].  Third,  because  it  main- 
tains a  feasible  design,  engineer  cannot  fail  to  meet  safety 
requirements  as  defined  by  the  contraints.  However,  the 
method  does  have  several  disadvantages  in  that  it  is  prone 
to  "zig-zag"  between  constraint  boundaries  and  that  it  is 
usually  does  not  achieve  a  precise  optimum.  This  method 
solves  the  nonlinear  programming  problem  by  moving  from  a 
feasible  point  (can  be  initially  infeasible)  to  another 
feasible  point  with  an  improved  value  of  the  objective 
value . 

The  following  strategy  is  typical  of  feasible  direction 
method  :  Assuming  that  an  initial  feasible  point  X^  is 
known,  first  find  a  usable-feasible  direction  S.  The  algo- 
rithm for  this  is  sinilar  to  linear  programming  and  comple- 
mentary pivoting  algorithms.  Having  found  the  search 
direction,  a  move  is  made  in  this  direction  to  update  the  X 
vector  according  to  Eg  (1.4)  .  The  scalar  a*  is  found  by  a 
one-dimensional  search  to  reduce  the  objective  function  as 
much   as   possible   subject  to   constraints.    That   is   IlIN 
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Figure  2.1    Algorithm  for  the  Feasible  Direction  Hethod. 

F(X  +  a*S)   subject  to  G(X  +  a*S)  <,    0.    It  is  assumed  that  the 
initial  design  xo   is  feasible,   but  if  it  is  not,   a  search 
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direction  is  found  which  will  direct  the  design  to  the 
feasible  region.  After  updating  the  xo  vector,  the  conver- 
gence test  must  be  performed  in  the  iterative  algorithm.  A 
convergence  criteria  used  in  this  is  implementation  are 
described  in  section  C.  The  general  algorithm  used  in  mscop 
is  given  in  Figure  2.1 

B.   SEARCH  DIRECTION 

In  the  feasible  direction  algorithm,  a  usable  -  feasible 
search  direction  S  is  found  which  will  reduce  the  objec- 
tive function  without  violating  any  constraints  for  some 
finite  move.  It  is  assumed  that  at  any  point  in  the  design 
space  (at  any  X)  the  value  of  the  objective  and  constraint 
functions  as  well  as  the  gradients  of  these  functions  with 
respect  to  the  design  variables  can  be  calculated.  Since 
these  gradients  cannot  usually  be  calculated  analytically, 
the  finite  difference  method  Eg  (2.1)  is  used  in  MSCOF. 

aF(X)      r  (X+  £e  )  -  F  (X) 

=  ~ (2.1) 


3X 
i 


where  e   is  the  ith  unit  vector 
i 

£  is  a  small  scalar. 

In  HSCOP,  £   is  0.1%   of  the  ith  design  variable 

In  the  feasible  direction  algorithm,  there  are  usually 
one  or  more  "active"  constraints.  A  constraint  G (X)  <  0  is 
"active"  at  X  if  g  (X)  ;:i  0.  As  shown  in  Figure  2.1,  if  no 
constraints  are  active  the  standard  steepest  descent  direc- 
tion  S  =  -  VF  is  used. 
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1.   Usable-Feasible  Direction 


; 

'2 

/ 
1 

""Tn,        ^Feasible    n. 
j^j/v^        sector            \^F(>0  =  const 

f\ 

> 

Usable-^ 
sector 

"""^^Ji/V^  Usable-feasible              \ 
^'>n:A.  secior            \                 \ 

\ 

'  -^1 

,  _.. 

Figure  2.2   Osable-Feasible  Direction 


Assume  there  are  NAC  active  constraints  at  X*   The  direction 
S,  is  "usable"  if  it  reduces  the  objective  function,  i.e.. 


VF«S  <  0 


(2.2) 


Similarly  the  direction  is  feasible   if  for  a  small  movement 
in  this  direction,  no  constraint  will  be  violated,  i.e.. 


7G • S  <  0 

This  is  shown  geometrically  in  Figure  2.2 


(2.3) 
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2 .  Active  Constraints 

It  is  necessary  to  determine  if  a  constraint  is 
active  or  violated  in  the  feasible  direction  algorithm.  A 
constraint  G  (X)  <  0  is  "active"  at  xo  if  G(XO)^0.  In  crder 
to  avoid  the  zigzagging  effect  between  one  or  more 
constraint  boundaries,  a  tolerance  band  about  zero  is  used 
for  determining  whether  or  not  a  constraint  is  active.  From 
the  engineering  point  of  view,  a  constraint  G (X)  <  0  is 
activ€  near  the  boundary  G(X)  =0  whenever  ACC  <  G (X)  <  VCC. 
ACC  is  the  active  constraint  criterion  and  VCC  is  the 
violated  constraint  criterion  in  .1SC0P.  Assuming  the 
feasible  constraints  are  normalized  so  that  G  (X)  ranges 
between  -1  and  0  for  reasonable  values  of  X,  the  constraint 
G(X)  <  0  is  considered  active  if  G (X)  >  -0.1.  The 
constraint  is  considered  to  be  violated  if  G (X )  >  0.004. 
This  is  an  algorithmic  trick  which  improves  efficiency  and 
reliability  of  the  algorithm.  However,  since  in  the  one  - 
dimensional  search,  all  interpolations  for  constraint  G  (X) 
are  done  for  zeros  of  a  linear  or  quadratic  appr ox imaticn  to 
G  (X)  in  crder  to  find  a*,  at  the  optimum  the  value  of  active 
constraints  are  very  near  zero,  but  may  be  as  large  as  0.004 
[Ref.  6].  From  an  engineering  point  of  view,  a  0.4  % 
constraint  violation  is  considered  to  be  acceptable. 

3 .  Subopt imizat  icn  Problem  and  Push-Off  Factors 

Zoutendijk   [Eef.  8]   has   shown   that   a   usable 
feasible  direction  S  nay  be  found  as  follows  : 

Maximize     fi  (2.  4) 

Subject  to  ; 

2F  (X)  -S  +  ^  <  0  (2.  5) 
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^G  (X)-  S  +  6.^   <  0 
S  bouEded 


j  e  J 


(2.6) 
(2.7) 


Where  scalar  ^  is  a  measure  of  the  satisfaction  of  the 
usability  and  feasibility  requirements.  The  scalar  Bj  in 
Eq(2.6)  is  referred  to  as  the  "push-off"  factor  which  effec- 
tively  pushes  the   search  direction   away   from  the   active 


AS 


g(\)  =  0 


-^X^ 


Figure  2.3    Push-Off  Factor  and  Bounding  of  the  S-Vector. 

constraints.  In  Eg  (2. 6)  ,  if  the  push-off  factor  is  7ero, 
the  search  direction  is  tangent  to  the  active  constraints, 
and  if  it  is  infinite,  then  the  search  direction  is  tangent 
to   the   objective  function.    It   has  been  found    that  a 
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push-off   factor      is    defined      as    follows      gives    good      results 
[Eef.    5:    p. 167]    : 


e     = 

J 


1    - 


G.{X)        1 

_J 

ACC 


e« 


(2.3) 


where    0^  =    1    . 

To  avoid  an  unbounded  solution  when  seeking  a  usatle 
-  feasible  direction  it  is  necessary  to  impose  bounds  on  the 
search  direction  S.  Cne  method  of  imposing  bounds  on  ssarch 
direction  is  to  impose  bounds  on  the  components  of  S-vector 
cf    form    : 


-  1    <    S.  <    1 


(2.9) 


This  choice  of  bounding  the  S-vector  actually  biases  the 
search  direction.  This  is  undesirable  since  we  wish  to  use 
the      push-off    factors  as      our      means    of      controlling      the 

search  direction.  A  method  which  avoids  this  bias  in  search 
direction  is  the  circle  as  shown  Figure  2.3  .  The  norm  here 
is 


S-S   <    1 


(2.9.1) 


4 .      Simple    Simplex-like    Method    for    Sea rch    Direction 

Vanderplaats  [Ref,  5:  pp. 168-169]  provides  the 
matrix  formulation  which  solves  the  above  sub-optimization 
problem    by    using    the    Zoutendijk    method. 


Maximize 


p.  v 


Subject    to    ; 


A.y    <  0 


(2.10) 


(2.11) 
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y.y    <    1 


(2.12) 


>Ihere 


y  = 


P   = 


(2. 13) 


n 


s 


A    = 


vG^(x),    e^ 


(2.14) 


T    ^ 
VF     (X),        1 

and    where    j   is    the    numhier    of    active    constraints     (NAC) 

When  the  solution  to  Ig(2.10)  through  (2.12)  is  found,  3  may 
be  normalized  to  some  value  other  than  unity,  but  the  form 
of   the   normalization    is    the    same.  A    solution    to    the    above 

problem  may  be  obtained  by  solving  the  following  system 
derived    from   the    Kuhn-Tucker    conditions    for    that    problem    : 


[e     x) 


=      c 


(2.  15) 


u      >    0  V      >    0  u-v    =    0 

i    -  i    -  ~    ~ 


(2.16) 


Where 
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B    =    -A- a"  (2.  17) 

I   =    Identity    matrix  (2.18) 

c   =    -A-P  (2.  19) 

Above  system  can  be  solved  using  a  complimentary  pivot  algo- 
rithm. Choose  an  initial  basic  solution  to  Ec(2.15)  is  to 
he 

v    =   c,  u    =    0  (2. 20) 

where  v  is  the  set  of  basic  variables  and  u  is  the  set  of 
nonbasic  variables.  If  all  v;_  >  0,  Eg  (2.  16)  is  also  satis- 
fied and  problem  is  solved.  If  some  v^  <  0,  the  solution 
procedure  is  as  follcvs  : 

Let  Ej,;,  be  the  diagonal  element  of  the  i-th  nonbasic  vari- 
able. 

1.  Given  the  condition  that  some  c  is  less  then  zero, 
we  find  max  {ci/Bn)  which  is  the  incoming  row  to  the 
basis. 
.  2.  The  incoming  column  is  changed  to  a  basic 
column,  the  tableau  is  updated  by  a  standard  simplex 
pivot  on  B[i   . 

3.  Until  all  c^>  0,  repeat  steps  1.  and  2. 

4.  When  all   c^  >  0,    the  iteration  is   complete.    The 
value  of  u  is  now  the  desired  solution. 

5.  3y  using     y  =   p-A -u,    we  get   the  usable-feasible 
search  direction   S   which  is  first  NDV  components  of 

y- 

^ •   Ifiiti^iil  l2^€asible  Designs 

The  method  of  feasible  directions  assumes  that  we 
begin  with  a  feasible  design  and  feasibility  is  maintained 
throughout  the  optimization  process.    If  the  initial  design 
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is  infeasitle,  then  a  search  direction  pointing  toward  the 
feasible  region  can  be  found  by  a  siniijle  modification  to 
direction  finding  problem. 

A  design  situation  can  exist  in  which  the  violated 
constraints  are  strongly  dependent  on  part  of  the  design 
variables,  while  the  objective  function  is  primarily  depen- 
dent on  the  other  design  variables.  This  suggests  a  method 
for  finding  a  search  direction  which  will*  simultaneously 
minimize  the  objective  while  overcoming  the  constraint 
violations.  These  considerations  lead  to  the  following 
statement  of  the  direction  finding  problem  [Ref.  5  : 
pp. 171-172]  : 


Maximize 


"  vr  (X)  •  S  +  5^ 


(2.21) 


Subject  to  ; 


vG  (X)  •  s  +  e  -^  <  0 

"     J    - 


je  J 


(2. 22) 


S-S  <  1 


(2.  23) 


where  J  is  the  set  of  active  and  violated  constraints,  and 
where  the  scalar  $  in  Eg (2. 21)  is  a  weighting  factor  deter- 
mining the  relative  importance  of  the  objective  and  the 
constraints.  Usually  a  value  of  5.  >  10000  will  ensure  that 
the  resulting  S-vector  will  point  toward  the  feasible 
region.  Incorporating  Eg  (2. 21)  and  Eg  (2.22)  into  the  direc- 
tion finding  algorithm  requires  only  that  we  modify  the 
p-vector  given  in  Eg  (2.  24)  and  the  A-matrix  of  Eg (2. 25). 


P  = 


VP(X) 

5 


(2.  24) 
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A  = 


~  3 


^i 


(2.25) 


0   <  50 


(2.26) 


Vie  use  the  simple  simplex-like  method  to  find  the 
search  direction  toward  the  feasible  region. 

C.   CNE-DIBENSIONAL  SEARCH 

T  •   No  Violated  Ccnstraints 

If  no  constraints  are  violated,  we  find  the  largest 
a*  in  Eci(1.U)  from  all  possible  values  that  will  minimize 
the  objective  on  S  without  violating  any  constraints,  active 
or  inactive. 

The  procedure  in  MSCOF  is  as  follows  : 

1.  Let  aO,  a1,  a2,  a3  be  the  scalar  in  11g(1.4)  corre- 
sponding to  pcints  xp ,  XJ,  X2/  XJ,  X4. 

2.  aC  =  0  at  given  point  XO. 

3.  In  order  to  get  a1,  we  can  calculate  the  a1  to 
reduce  the  objective  by  at  most  10%  or  to  change  each 
of  the  design  variable  X  ^Y    ^t  most  ^Q%. 

4.  Update  the  design  variables  to  XJ  '^sing  Eg  (1.4)  . 

5.  Evaluate  the  objective  for  XJ,  and  check  the  feasi- 
bility. If  one  or  more  constraints  is  violated,  then 
a1  is  reduced  to  a1/2,  and  we  go  to  step  4. 

6.  In  order  to  estimate  a2,  we  can  use  the  quadratic 
approximation  with  2  points  X,    XJ  and  the  VE. 
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7.  Update  the  design  variables  to  X2  by  Eq(1.4)  and 
check  the  side  constraints. 

8.  Evaluate  the  objective  and  constraints. 

9.  Now  having  3  a's,  and  values  of  objectives  and 
constraints  for  design  variables  XO,  XJ ,  X2  are 
known,  so  by  using  3-point  quadratic  approximation,  a 
value  of  aS  is  found. 

10.  Update  the  new  optimal  point  in  search  direction  by 
Eg{1.U)  . 

11.  Evaluate  the  objective  and  constraints. 

12.  Now  choose  last  3  values,  a1 ,  a2,  a3  and  find  a  new 
a3  using  3-points  Quadratic  approximation 

13.  Choose  the  a*  among  the  5  points  which  corresponds  to 
the  minimum  objective  function  value  with  no-viclated 
constraints. 

2  .   On e  or  More  Const raints  Viola ted 

If  one  or  more  constraints  are  initially  violated,  a 
modified  usable-feasible  direction  is  found.  It  is  then 
necessary  to  find  the  scalar  a*  in  Eg  (1.4)  which  will  mini- 
mize the  maximum  constraint  violation,  using  the  most 
violated  constraint  j,  a  good  initial  estimate  for  a*  is 

-G  .  (X) 

a*  = (2.27) 

VG  (X)  .  S 
3 

Since  the  gradients  of  the  violated  constraints  are 
known,  the  scalar  which  is  required  to  obtain  a  feasible 
design  with  respect  to  violated  constraint  in  the  search 
direction,  is  given  to  a  first  approximation  by  Eq(2.27). 

The  more  detail  procedure  in  :iSCOP  is  as  follow  ; 

1.  Choose  the  most  violated  constraint  j. 

2.  Calculate  a*  for  violated  constraint  j  using 
Eg  (2.27)  . 
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3.  Update  the  design  variables  for  a*  and  check  the  side 
constraints. 

4.  If  one  or  more  violated  constraints  still  exist,  then 
calculate  the  derivative  of  objective,  violated  and 
active  constraints  and  find  a  new  search  direction 
and  then  go  to  step  1.  Otherwise  proceed  with  the 
optimization  in  the  normal  fashion. 

E.   CCNVEBGENCE  CRITIBIA 

A  desired  property  of  an  algorithm  for  solving  a  nonli- 
near problem  is  that  it  should  generate  a  sequence  of  points 
converging  to  a  global  optimal  point.  In  many  cases, 
however,  we  may  have  to  be  satisfied  with  less  faverable 
outcomes.  In  fact,  as  a  result  of  non-convexity,  problem 
size,  and  other  difficulties,  we  may  stop  the  iterative 
procedure  if  a  point  belongs  to  a  described  set,  which  is 
defined  in  HSCOP  as  fellows  ; 

1.  Q  =  1 X  1  IK°  -  II  <  £xlxo|  } 

2.  Q   =  {X     \     !F(XO)  -  F{X)1  <   e-|F(XO)|  } 

2  ^ 

In  MSCOF,  the  algorithm  is  terminated  if  a  point  _X  is 
reached  such  that  X  ^  Q,  f)  Q^  -  £,.  is  0.001  and  Sp  is 
approximatly  0.001.  Since  in  engineering  design  problems  it 
is  not  necessary  to  find  solutions  with  more  than  three 
significant  digits. 
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III.     HSCOP    USAGE 

A.  INTEODDCTION 

Since  this  MSCOP  is  written  in  WATERLOO  BASIC  Version 
2.0,  it  is  verv  convenient  to  use.  The  user  must  first 
formulate  the  design  problem  with  the  classical  machine 
design      criteria.  Given      the    formulation      of      the      design 

problem  as  a  nonlinear  program,  the  user  then  enters  the 
problem  as  a  part  of  a  BASIC  program.  The  user  defines  the 
objective  function  and  constraint  functions  using  EASIC 
statements.  Other  parameters  are  input  as  data  :  the  number 
of      design        variables      NDV,  the      number        of      inequality 

constraints  NIQC,  variable  bounds  an  initial  design  value 
and    a    print   control    number. 

B.  PECBIEM    FOBMDLATICN 

Generally,  the  experienced  design  engineer  will  be  able 
to  choose  the  appropriate  objective  for  optimization 
depending  on  the  requirements  of  the  particular  application. 
The  physical  phenomena  of  significance  should  first  be 
summarized  for  the  device  to  be  designed.  The  appropriate 
objective  can  then  be  selected  and  constraints  can  be 
imposed  en  the  remaining  phenomena  to  assure  an  acceptable 
design  from  all  standpoints.  However,  the  initial  formula- 
tion for  the  optimization  problem  should  not  be  more  compli- 
cated then  necessary  and  this  often  requires  the  making  of 
some  simplifying  assumptions.   [Eef.  9]. 

After  completing  the  formulation  of  the  design  prcblem, 
the  design  engineer  should  be  able  to  answer  the  following 
questions  : 

1.   What  are  the  design  variables  ? 
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2.  What  is  the  objective  function  ? 

3.  What  are  the  inequality  constraints  ? 

4.  What  are  the  hounds  on  the  variables  ? 

The  engineer  is  then  ready  to  input  the  program  to  the 
MSCOF.  However,  additional  study  and  preparation  of  the 
problem  may  be  useful.  In  particular,  redundant  constraints 
should  be  avoided  if  possible.  MSCOP  will  operate  with 
redundant  constraints  but  it  will  operate  faster  without 
them.  Selection  of  an  initial  design  point  from  which  to 
start  this  program  is  important,  since  it  affects  perform- 
ance and  running  time.  The  user  should  use  any  available 
information  which  gives  a  good  initial  approximation.  If 
side  constraints  exist,  the  user  must  be  sure  the  initial 
values  of  the  design  variables  do  not  violate  the  side 
constraints.  This  program  will  automatically  handle  an 
initial  design  point  which  is  infeasible  with  respect  to  the 
G(X)  <  0  constraints.  However,  if  the  initial  point  does  not 
violate  these  constraints,  convergence  will  likely  be  more 
rapid . 

C.   PBCBIEM  ENTEY 

Problem  entry  is  accomplished  by  editing  the  main 
program  directly.  As  an  example,  consider  the  following 
simple  NIP  with  two  design  variables,  and  three  constraint 
functions. 

2  2 

Minimize    F(X)  =  X  +  3XX   +2X-X-X+1 

1      12       2    12 


subject  to  ; 


X   +  X   -  3  <  0 
1     2 


1      1 

—  + 2  <  0 

XX- 
1      2 
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2 
X   +X-X-2<0 
1     12- 


X   >  0.  1 
i  - 


With  the  MSCOP  loaded  into  memory  and  listed  on  the  CRT, 
modifications  are  made  on  the  program  lines  as  follows  to 
input  this  example  : 

line  100 

Just  after  the  word  "data",  three  integers  are  added, 
separated  by  a  conma.  The  first  number  is  ND V  vhich  is 
the  number  of  design  variables,  the  second  is  NIQC  which 
is  the  number  of  ireguality  constraints,  and  the  third  is 
IPET  which  is  print  control  number  (  0  ;  only  final 
results,  1  ;  given  data  and  final  results,  2  ;  given  data 
and  iterative  subcptimal  results) 

for  example  : 

100  data  2,3,2 

Lines    201-220 

Each  line  here  corresponds  to  a  separate  design  variable, 
beginning  with  X  (1)  and  continuing  in  order  to  input 
X  (NDV) .  On  each  line,  three  values  are  separated  by 
commas.  After  the  word  "data",  these  values  are  the 
initial  values  of  the  design  variable,  the  lower  bound  on 
the  variable  and  the  upper  bound  on  the  variable.  If  no 
bound   is   to    be   specified,    the    entry   is    filled    by    "no". 

For   the   sample   problem,    the   input   is    : 

201  data    3.  ,0.  1,no 

202  data    3,, 0.1, no 
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lines  400  -  450 

These  lines  are  available  for  defining  the  objective 
function.  The  objective  function  must  be  defined  in 
terms  of  subscripted  design  variables  X  (  1) ,  X {2)  ,    etc. 

Tor  the  sample  protlem,  the  input  is  : 

400  fn_f  =  X  (1)**2  +  x(1)*x(2) +2.*x  (2)  **2-x  (  1) -X  (2)  +1. 

lines  500-650 

These  lines  are  available  for  defining  the  inequality 
constraint  functions,  which  must  be  expressed  using  the 
format  : 

601  if  i  =  k  then  fn_g  =  G  (x) -  b 

i      i 

For  the  sample  problem,  the  input  is  : 

OC601  if  i  =  1  then  fn  g  =  x(1)+x(2)-3. 
00602  if  i  =  2  then  fn  g  =  1 . /x  ( 1)  +  1  ./x  (2) -2  . 
G0603  if  i  =  3  then  fn^g  =  x  ( 1)  **2  +  x  (1 ) -x  (2) -2.  - 

If  there  are  many  constant  values  in  the  constraint  func- 
tions, the  user  may  input  data  for  these  functions  on 
lines  501-600  in  order  to  simplify  their  statements. 
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IV.  EXAHPIE  PROBLEMS 


A.   DESIGN  OF  CAHTILE7ERED  BEAM 


1.   Uniform  Cantilevered  Beam 


Assume  a   cantilevered  beam  as   shown  in   Figure  4. 1 
must  be   designed.    The   objective  is   to  find   the  minimum 


1  =  20O"' 


r 


T 

H 


(J 

=  20000  Psi 

E 

=  30  E  6  Psi 

y 

=1.0  inch 

p 

=  10000  lbs 

Figure  4.1    Design  of  a  Dniforn  Cantilevered  Beam. 


volume  of  material  which  will  support  the  load  P. 

The  design  variables  are  the  width  B  and  height  H  in 
the  team.  The  design  task  is  as  follows  :  Find  B  and  H  to 


irinimize  volume    V  =  B  H  1 


(U.I) 
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ve  wish  to  design  the  beam  subject  to  limit  on  bending 
stress,  shear  stress,  deflection  and  geometric  conditior.s. 
The  bending  stress  in  the  beam  must  not  exceed  20,000  psi . 

«  c      6  P  1 

Cr  =  =  <  20,000  (4.2) 

b      I  2    - 

B  H 

The   shear    stress    must   not   exceed    10,000    psi. 

^  3    P  3    P 

(X=    =    <    10,000  (U.S) 

•"  2    A  2    B    H      -  ' 

and   the    deflection    under    the    load   must   not    exceed    1    inch. 

3                            3 
C          PI                  4    P   1 
O   =    =    <    1.0  (4.4) 

3    E    I  3         - 

E    B    H 

Additionally,  geometric  limits  are  imposed  on  the  beam  size. 
0.5<B<5.0  (4.5) 


1.0  <  H  <  20.0  (4.6) 

H/b  <  10.  (4.  7) 

Now  we  can  input  this  problem  to  MSCOP. 

Input  NDV,  NIQC,  IPRT 

00100  data  2,4,2 

Initial  starting  points 

00210  data  3.5,0.5,5.0 
00220  data  16.0,1.0,20.0 

Evaluation  of  objective 

00400  fn_f  =  tl*x(1)*x(2) 
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Evaluation  of  constraints 
00500  tl  =  200. 

30.e+6 
10000. 

1  then  fn_g 


00SG1  hG 

00S0  2  bp 

00503  if 

00S03  if 

00503  if 

00503  if 


then  fn'_g 
then  fn_g 
then  fn_g 


=  6.  *bp*tl/f20000  .*b*h*''  2) -1 
=  3. *bp/( 10000.*2.*b*h) -1. 
=  U. *bp*tl**3/ (be*b*h**3) -1. 
=  h/b-10. 


TABLE  I 
The  Solution  of  a   Uniform  Cantilevered  Beam 


objective   ;   6664.0 

design  variable  ; 

X(1)  =  1.852 
X{2)  =  17.99 

constraint  ; 


9(2) 
g(3) 

gC*) 


0.000902 
-0-9549 
-0.0109 
-0.0286 


As   a    result   of    this    problem    are    in    Table    U.I. 

2  .      ^Variable  Can  tilever ed    Beam 

The  cantilevered  beam  shown  in  Figure  4.2  is  to  be 
designed  for  minimum  material  volume.  The  design  variables 
are  the  width  b  and  height  h  at  each  of  5  segments.  We 
wish  to  design  the  beam  subject  to  limits  on 
stress  (calculated  at  left  end  of  each  segment),  deflection 
under  the  load,  and  the  geometric  requirement  that  the 
hoight    of   any    segment   does    not    exceed   20    times    the    width. 
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Cross  section 


Figure  U.2    Design  of  a  Variable  Cantilevered  Beam. 

The  deflection  y    at  the  right  end  of   segment  i  is 
calculated  by  the  following  recursion  formulas  : 


y   =  y   =  0 
o     o 


(4.8) 


P  1 


I   I 


1     1 

L  +  —  +  s:  1 

2     j=1   i 


i-1 


{a.9) 


p  1 


y  = 


2  E  I 


2  1 
i  i 

L  -  ZT  1   +  

j=1   i     3 


+  y.  1.  +  y.  ,   (4.10) 
1-1  1   1-1 
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vhere    the    deflection    y   is      defined    as   positive    downward,  y' 

is    th€    derivative    of    y    with    respect    to    the    X,      and    1;.    is  the 

length    of   of      segment  i.      Young's    raodjlus    E    is      the    same  for 
all   segments,    and   the   moment   of   inertia   for    segment    i    is 


I      = 
i 


3 

t      h 
i      i 

12 


(U.  11) 


The  tending   moment  at  the  left   end  of  segment  i   is  calcu- 
lated as 


1   =PL  +  1   -H   1 
i      I      i    3=1   i  J 


(4.12) 


and  the  corresponding  maximum  bending  stress  is 


0"  = 

i 


M   h 
i   i 

2  I 


('4.13) 


The  design  task  is  now  defined  as 


Minimize 


Subject  to  : 


V  =  z:  b  h  1 
i=1   i   i   i 


(U.  14) 
(4.15) 


(Ti 


-  1  <  0 


i  =  1,  .  .  . , N 


(4.  16) 


N 


-  1  <   0 


(4.  17) 


h   -  20  b   <   0 
i       i  - 


i  =  1,...,N 


(4.  18) 
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b   >  1.0 
i  - 


h   >  5.0 
i  - 


i  =  1,  .  .  .,N 


(U.19) 


where  ?  is  the  allovable  bending  stress  and  y  is  the  allo- 
wable displacement.  This  is  a  design  problem  in  10  vari- 
ables. There  are  6  nonlinear  constraints  defined  by  Eq  (4. 16) 
and  Eg  (4. 17),  and  5  linear  constraints  defined  by  Eg  (4. 18), 
and  10  side  constraints  on  the  design  variables  defined  by 
Eg(4.  19)  . 

Now  we  can  input  this  problem  to  ^SCOP. 

Input  NDV,  NIQC,  IPRT 

00100  data  10,  1  1,2 

Initial  starting  points 

0021 0  data  5. , 1 .  ,no 
00220  data  5. , 1 .,no  - 
00230  data  5. , 1 . ,no 
00240  data  5.,1.,no 
00250  data  5. , 1. ,no 
00260  data  40. ,5. , no 
00270  data  40. ,5. , no 
00280  3ata  40. , 5. ,no 
00290  data  40. ,5. , no 
00300  data  40. ,5. ,no 


Evaluation  of  objective 

00400  fn_f  =  100.   *  (   x(1)*x(6) 
x(4)  *x(9)  +  x(5)  *x  (10)  ) 


+  X  (2)  *x  (7)   +  X  (3)  *x  (3) 


Evaluatiom  of  constraints. 


00490  def  fn  q  {x,i\ 
00500  pcb  =^^0006 


00498  dim  bmT 


bi  (10)  ,sigi  ( 10)  ,y  pb  (10)  ,  yb  (10) 


00501  be  =  200. e+5 

00502  tl  =  200. 

00503  sigb  =  14000. 
005C4  ytb  =  .5 

00505  si  =  40. 

00506  for  m  =  1  to  5 

00507  bm(m)  =  pcb*  (tl+sl-m*sl) 

00508  next  m 

00509  for  ffi  =  1  to  5 

00510  km  =  m+5 

00511  bi(m)  =  x  (m)  *x  (km)  **3/1  2. 

C0512    sigi (m)  =  bm (m  *x (km) /  (2. *bi  (ra) ) 

00513  next  m 

00514  yzo  =  0. 

00515  yp2  0  =  0. 

00516  for  m  =  1  to  5 
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00517 
00S18 
00519 
0  0  5  2  0 
00521 
00522 
00550 
00560 
00570 
00560 
00590 
00600 
00610 
00620 
00630 
00632 
0063a 
00636 


y 
y 

next 

rem 

if  i 
i 
i 
i 
i 
i 
i 
i 
i 
i 
i 


pb  (m 
umm 
L  (ai) 
pzo 
zo  = 


)  =  (pcb*sl* (tl+sl/2.-m*sl) ) /  (be*bi  (m) ) -ypzo 
=  pcb*£l**2* (tlrra*sl+2.*sl/3.) 


=  dumm/(2 
=  ypb  (E) 
yb(ra) 


*be*ti  (m) )  +ypzo*sl  +  yzo 


if 
if 
if 
if 
if 
if 
if 
if 
if 
if 


m 
constraint  function 
1  then  fn_g  =  si 
si 
si 
si 
si 


1  then  fn_g  = 

2  then  fn_g  = 

3  then  fn_g  = 

4  then  fn_g  = 

5  then  fn_g  = 

6  then  fn_g  =  yb 

7  then  fn_g  =  x  ' 

8  then  fn_g  =  x 

9  then  fn_g  =  x 

10  then  fn_g  =  x 

11  then  fn_g  =  x 


sigb- 

sigb- 

sigb- 

sigb- 

siqb- 
b-1. 

*x  (1) 

*x  2 

*x  3 

.*x  ( 
0.*x 

\h 

TABLE  II 
The  Solution  of  a  Variable  Cantilevered  Bea: 


objective   ;   62133-35 

design  variables 

X(1)  =  2.994 

X(2)  =  2.782 

X(3)  =  2.528 

X(4)  =  2-203 

X(5)  =  1.761 

X(6)  =  59.88 

X(7)  =  55.62 

X(8)  =  50.56 

X(9)  =  44.14 

XOO)  =  35-  19 


constraints 


G(1)  = 
G(2)  = 
G(3)  = 
G(4)  = 
G(5)  = 
G(6)  = 
G(7)  = 
G(8)  - 
G{9)  = 
G{10)  = 
G(11)  = 


-0.00219 
-0.00415 
-0-00508 
-0-00406 
-0-0177 
-0- 4401 
-0.0101 
-0.0231 
0.0000 
-0.0248 
-0-0278 
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E,   SIHPIE  TRUSS 


Y 

4 

< 

P  =  50000  N 

E  =  200  GPa 

(J  =  +  14000  N/cm 

t  =  100  cm 

^x^ 

4 

I 

1 

© 

< 

/ 

t 

^ 

3 

1    - 

2     '^ 

/ 

^ 

P 

Figure  U. 3   Design  of  a  5-Bar  Truss. 

A  simple  truss  with  5  members  as  shown  in  Figure  U.3  is 
designed  for  the  minimum  volume.  The  design  variables  are 
the  sectioEal  areas  of  the  members.  The  constraints  are 
formed  for  the  stresses  of  the  members  not  to  exceed  the 
given  allowable  stress.  The  lower  bound  for  each  design 
variable  is  also  considered.  The  stresses  are  obtained  by 
the  displacement  method  of  the  finite  element  analysis.  The 
equation  to  be  solved  is  given  by 


K-u  =  P 


(4.20) 


where   K  is   the  stiffness   matrix,   u   is  the   displacement 
vector  and  P  is  the  lead  vector  as  follows  : 
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u  = 


r                   \ 

u 

1 

V 

1 

u 
2 

V 

2 

p  = 


0 

0 

0 

-5000 


(4.21) 


K  =E 


A.             ^5 

As 

0 

0 

A5 

A,      P^s 
X"  V5fi. 

0 

/la 

0 

0 

A, 
JL 

0 

- 

A4. 

A2.      f^'^ 

(4.12) 


From  Eq.  (4.20)  the  displacements  are  solved  by 

-1 

D  =  K  -P 


(4.23) 


Having    displacements      at    all      nodes,      we      can   calculate      the 
stress   for   each    element. 


E'Ai 


01=  E-g  = 


(4.24) 


where 


/  2  2 

Al      =      1(1    +   u    )       +    V 
1^11  1 


2  2 

Al      =      /(l+v-v)    +    (u-u)-l 
2  V        2         1         2  1         2  2 


=  / 


2  2 

Al=/(l+u)        +v        -       1 
3        v'       3        2  2  3 


(4.25) 
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^1       =    /(     1     +     U     )         +      {     1    -     V    )         -     1 

U        J         3        2  2         2  a 


2  2 

Al      =/(l+u)        +(l    +    v)       -1 
5        >i        3         1  2         1  5 


The    design    problem    is  given    by 


minimize  V   =     21      A      1 

i=1      i      i 


Subject    tc 


-1.0<0      i=1,...,5 


(4.26) 


(4.27) 


A       >    0.  1 
i    - 


i    =    1 , . .  .  ,  5 


(4.28) 


The    MSCOF   input    for    this    problem    is    given    as    follows 
Input    NDV,    NIQC,    IPRT 

00100    data    5,5,2 
Initial    starting    point 


00200  data  3. , . 1,no 
00202  data  3. , .  1,no 
00204  data  3. , .  1,no 
00206  data  3. , . 1,no 
00208    data    3. , .  1,no 


Evaluation    of  objective 

00400    fn    f    =  100    *     (   x(1)     +    X  (2)       +    x(3)       +    sqr(2.)*x(4) 
sqr  (2.)  *x75)  ) 

Evaluation    of  constraints 


0500    dim    vv  (5) 
050  1    te   =    2.e+7 

0502  tl  =  100. 

0503  sigb  =  14000. 
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0504 
050S 
0506 
0507 
0506 
0509 
0510 
051  1 
0512 
0513 
0514 
0515 
0516 
0517 
0518 
0519 
0520 
0  52  1 
0522 
0523 
0590 
0592 
0594 
0595 
0596 
0598 
0600 
0602 
0604 
0606 
0608 
0610 
06  50 


sqr  (2.) 

1) +x  (5)/cs)  *ct 

) *Ct/C£ 

2)+x  {5)/cs)  *ct 
2)  *ct 

3  +x  (4)/cs)  *ct 
4) *ct/cs 


4)/cs)  *ct 

2*dk1)  /lc24 
/k34 


f  nend 


raint 

t  hen 

then 

then 

then 

then 


2*dk1+k4  3*ak3+k4  4*dk2) 

K 

vv  (  1 

vv(2 

vvh 

1 

vv 

vv 


) **2  +  vv  (2) **2) -tl 
-vvf4))  **2+  (vv(11-vv(3) 
)**2  +  vv  (4)  **2)  -tl 

(3)  )  **2+  (hl-vv 
(1)  )**2+  (hl  +  vv 


)  ♦*2)  -tl 


(4)  )**2)  -hi 
(2)  i**2i-hl 


f  n_g 
f  n_g 
f  n_g 

f  n_g 
f  n  g 


te*dl1/ 

te*dl2/ 
te*dl3/ 
te*dl4/ 
te*dl5/ 


tl*sigb 

tl*sigb 

tl*sigb 

hl*sigb) 

hl*sigbi- 


TABLE  III 
The  Sclution  of  a  5-Bar  Truss 


objective   ;   108-52 

design  variables 

X(1)  =  0-  1 

X (2)  =  0.  1 

X(3)  =  3.514 

X(4)  =  4.948 

X(5)  =  0.1 


constcaint 


G(1)  = 

G(2)  = 

G{3)  = 

G(4)  = 

G(5)  = 


-1.9988 
-2.0030 
-0.0030 
-0.  1203 
-1.8797 
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V.  SDH MARY  AND  CONCLnSION 

Numerical  optimization  is  a  powerful  technique  for  those 
confronted  with  practical  engineering  design  problems.  It 
is  also,  a  useful  tool  for  obtaining  reasonable  solutions  to 
the  classical  engineering  design  problems.  Since  many  engi- 
neers are  now  using  nicrocomputers  for  solving  design  prob- 
lems, the  development  of  microcomputer  software  which  can  be 
easily  used  is  needed. 

In  this  thesis,  an  algorithm  for  constrained  optiiriza- 
tion  problems  is  programmed  in  standard  BASIC  language 
(WBASIC  version  2.0)  on  an  IBM  3033.  The  users  can  easily 
convert  this  to  other  nicrocomputers. 

MSCOI  (Microcomputer  Software  for  Constrained 
Optimization  Problems)  employs  the  method  of  feasible  direc- 
tions and  specific  modifications  of  a  one-dimensional  search 
for  constrained  optimization.  MSCOP  has  been  validated  by 
tests  on  three  constrained  optimization  problems.  Its 
performance  is  good  and  could  be  made  better  through  refine- 
ment cf  the  algorithm. 

Since  microcomputers  are  available  with  reasonable 
memory  size  and  computational  speed,  their  capabilities  will 
continue  to  improve  as  more  engineering  software  becomes 
available.  MSCOP  is  considered  to  be  a  first  step  toward 
more  widespread  use  cf  optimization  techniques  on  microccm- 
puters. 
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APPENDIX  A 
MSCOP  PHOGEAM  LISTING 


0010 
0020 
00  21 
0030 
0040 
0050 
0060 
0070 
0030 
0090 
0100 
0110 
0115 
0120 
0125 
0130 
0135 
0140 

0150 

0160 
0200 
0210 
0360 
0370 
0375 
0380 
0390 
0400 
0410 
0420 
0430 
0440 
0450 
0480 
0490 
0500 
0510 
0520 
053  0 

0540 
0550 

0560 
0650 
0700 
0710 
0720 
0730 
0740 
0750 
0760 
0770 
0780 
0785 
0790 
0800 


cptio 

dllE     X 

din! 
dim 
dim 
dim 
rem 
gosufc 
rem  i 
read 
data 
for  i 
rem 
re 
xO 
if 
re 
if 
va 
if 
va 
next 
data 
data 
rem  e 
ot  j  = 
itri 
rem  o 
def  f 
fn 
f  nend 
rem  e 
for  i 

next 
rem  c 
def  f 
tl  = 
be  = 
bp  = 
2.1  i 


n  bas 

(21)  L 

Hta  (^ 
(51,2 
wrk  (5 
rku(5 
npu  t 

1000 
nput 
ndv ,  n 
2,4^2 

input 
ad   x  ( 

(il    = 

niqc 
ad    lo 

lo?; 

lue  (1 

ups 
lue  (u 

3.5,0 
16.  ,  1 
valut 

fn    f 
=    1" 
bject 
n   f  (X 

f   = 


e    1 
xO  (21 


(21),qcv(51J  ,ngcv{51)  ,  df  (21)  ,  dg  (51 ,  2  1) 
,  wrxY  (51,  5  1)^ 

lj  ;  wrki  (5lJ  ;iowb  (21)  ;uprb  (2  1)  ;ioS  {6)  ,up?  (6) 


1)  ;wrkyt5i;51 

1    ,b(5l,5lf  ,p(21),y(21),sf21|,u( 

1) , :vrk    51) , wrki  (51)  ,wrK2  (51) ,wr 


i^ 


51)  ,c  (51) 
k3(§1) 


data 

0 

number  of   design    variables    and    constraints. 

iqc, ifrt 

to  n  dv 
initial  value  of  design  variables 

i)     . 
x(i) 

=    0    then    160 
S^upl 
-    *no'    then    lowb  (i)     =    bnlo    else    lowb(i)    = 


o$) 

=  'no'  then  uprb  (i)  =  bnup  else  uorb  (i)  = 

p$) 

.5,10. 
.0,20. 
e  the  objective-function 

(X) 

ive    function 
200.*x(1)*x  (2) 


te    the   constraints 

to    niac 

=    fn_g(x,i) 

aint    functions 
,i) 


if    i    =   2    t 
if    i   =   3    t 


valua 
=    1 

Y(i) 

1 

onstr 

n   g  (X 

2Ua. 

30.e+6 

10000. 

=    1    then    fn_g    = 


hen    fn_g    = 
hen    fn_g    = 


if  i 
f  nend 
rem  i 
ical 
if  ic 
rem  c 
gosub 
rem  p 
rem 
rem  r 
ical 
if  ic 
rem  1 
for 


=   4    t 


6.  *bp*tl)/ 

'2000O.*x  (1)  *x  (2)  **2)  -1  . 
'3.  *bp)/(200  0  0.*x  (1)  *x  (2)  )  -1. 
'4.*bp*tl**3)/ 
■be*x  M)  *x  (2)**3)  -1. 


hen    fn_g   =    x  (2)  /  ( 10  .  *x  ( 1)  )  -  1 . 
1  counting   number    input 


nitia 
=    1 
al    > 
all    t 
2000 
rint    results. 


3  then  stop 

he  optimization  code, 


e-counting  number  input. 

=  ical+1 

al  =  3  then  850 

0%   reduce  the  design  variables. 

i  =  1  to  ndv 
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0810     x(i)  =  0.9*x(i) 

0320    xO  (i)  =  X  (i) 

0830   next  i 

0840  goto  720 

0850  rem  lOf  increase  design  variables. 

0860  for  i  =  1  to  ndv 

0870    X  (i)  =  1.1*x(i) 

0880     xO  (1)  =  X  (i) 

0890  next  i 

0900  goto  720 

2000  rem  calculate  the  obi.  constraint  fen. 

2001  obj  =  fn  f  (X) 

2002  for  i  =  T  to  nigc 

2003  2cv(i)  =  fn_g{x,i) 

2004  next  i 
2003  itrq  =  1 
2010  itrq  =  itrg+1 

2020  rem  calculate  the  number  of  active  and  violate 

constraints. 
2030  gosub  3500 
2040  rem  calculate  the  gradient  of  objective  and 

active  or  violated  constraints. 
2050  gosub  3800 
2060  if  nave  =  0  then  2190 
2070  gosub  3900 

2080  rem  calculate  the  push-off  factors 
2090  gosub  UOOO 
2100  rem  making  the  latrix  c 
2110  rem  normalized  the  df (i) 
2120  gosub  4100 
2130  rem  normalized  the  DG (i) 
2140  gosub  4200 

2150  if  nvc  >  0  then  gosub  4400  else  gosub  4600 
2160  rem  calaulate  the  usable-feasible  direction  s  (i) 
2170  gosub  5000 
2180  goto  2230 
2190  rem  normalize  the  df(i) 
2200  for  i  =  1  to  ndv 
2210     s(i)  =  -(df  (i).) 
2220  next  i 

2230  rem  normalize  the  s (i) 
2240  gosub  5700 

2250  rem  one-dimensicnal  search 

2260  if  nvc  =  0  then  gosub  6000  else  gosub  9000 
2270  rem  update  x  for  alph 
2280  gosub  7000 
2290  gosub  7100 

2300  rem  calculate  new  point  value. 
2310  nobj  =  fn_f  (x) 
2320  rem  convergence  test 
2330  gosub  6780 

2340  if  walp  <=  accx  and  delf  <=  dabf  then  2470 
2350     itri  =  itri+1 

2360     if  itri  >  mxit  then  print  'check  the  problem' 
2370     obj  =  nobj 
2380  for  i  =  1  to  ndv 
2390     xO  (i)  =  X  (i) 
2400  next  i 

2410  for  i  =  1  to  nice 
2420     gcv{i)  =  fn  g(x,i) 
2430  next  i        -^ v  '  / 
2440  if  iprt  =  2  then  2460 
2450     gosub  9200 
2460  goto  2010 

2470  rem  print  final  results 
2480  print  '*****  final  results  *****  ' 
2490     gosub  9200 
2500  return 
3000  rem  initialize  the  integer  working  array 
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3005    for    i   =    1    to    niqm 

3010  iwrk(i)    =    0 

3015   next    i 

3020    return 

3050    rem   initialize    the    integer    working    array 

3055    for   i   =    1    to    nicjin 

3060  jwrk(i)     =    0 

3065    next    i 

3070    return 

310C    rem   initialize  the    one-dimension    working    arrav 

3105  for  i  =  1  to  niqm 

3110     wrki  (i)  =0. 

3115  next  i 

3120  return 

3150  rem  initialize  the  one-dimension  working  array 

3155  for  i  =  1  to  niqm 

3160     wrk2(i)  =  0. 

3165  next  i 

3170  return 

3200    rem    initialize    the    one-dimension   working    array 

320  5    for    i   =    1    to    niqc 

3210  wrk3  (i)    =    gcv(i) 

3215    next    i 

3220    return 

3250    rem   initialize    the    two-dimension    working    array 

3255    for   i   =    1    to    niqm 

3260  for    i    =    1    tc  ndvm 

3265  wrky(i,j)    =    0. 

3270  next    j 

327  5    next    i 

3280    return 

3300    rem   initialize    the    derivative    of   objective    DF  (i) 

3305    for    i   =    1    to    ndvm 

3310  df  (i)     =   0. 

3315    next    i 

3320    return 

3350    rem   initialize    the   a  (i  ,  j)  ,  p  (i)  ,y  (i)  ,c  (i) 

3353    for    i   =    1    to    ndvm 

3356     p  (i)  =  0. 

3359     y(i)  =  0. 

3362     for  j  =  1  tc  niqm 

3365       a(^,i)  =0. 

3368     next  j 

337  1  next  i 

3374  for  j  =  1  to  nicm 

3377    cij)    =    0. 

3380  next  j 

3383  return 

3400  rem  initialize  the  derivative  of  constraints  DG(i,j) 

3405  for  i  =  1  to  niqm 

3410    for  1  =  1  tc  ndvm 

341^       dg(i,j)  =  0. 

3420     next  3 

3425  next  i 

343  0  return 

3450  rem  initialize  the  b(i,j) 

3455  for  i  =  1  to  niqm 

3460     for  1=1  to^niqm 

34  6  5        b{i,j)  =0. 

3470     next  j 

3475  next  i 

3480  return 

3500  rem  Calculate  the  number  of  active  and  violate 

constraints. 
3502  gosub  3000 
3504  gosub  3100 
3510  nac  =  0 
3520  nvc  =  0 
3530  for  i  =  1  to  nigc 
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3540 
3550 
3560 
3570 
3530 
3590 
3610 
3620 
3630 
3640 
3650 
3660 
3670 
3680 
3690 
3700 
3710 
3720 
3730 
3740 
3750 
3790 
3800 
3805 
3810 
3815 
3820 
3825 
3830 

3  83  5 
3840 
3850 
3860 
3900 
3901 
3910 
3915 
3920 
3925 
393  0 
3935 
3940 
3945 
3950 
3955 
3960 
3966 
4000 
4010 
4020 
4030 
4040 
4090 
4100 

4  102 
4105 
41  10 
4115 
4120 
4125 
4127 
4130 
4135 
4140 
4145 
4200 
4202 
4205 
4210 


if  gcv  (i)  >=  vcc  then  3580 

if  gcv  (i)  <  ace  then  3590 
nac  =  nac+1 

goto  3590 

nvc  =  nvc+1 
next  i 

nave  =  nac+nvc 
if  nave  =  0  then  3790 

ii  =  1 


ij.=  1 

ror,i  =  1  tc  nigc 


if  gcv(i)  >=  "vcc  then  3720 
if  gcv  (i)  <  ace  then  3750 

iwrk.  (nvc+ii)  =  i 

wrki  (nvc+ii)  =  gcv  (i) 

ii  =  ii+1 
goto  3750 
iwrk  (j j)  =  i 
wrki  (^ j)  =  gcv (i) 

Jl=  3  3+ 1 

next  1 
return 

rem  calculate  the  gradient  of  f  (x) 
gosub  3300 
for  i  =  1  to  ndv 
dxi  =  fdm*abs  (X  (i)  ) 

if  dxi  <=  mfds  then  dxi  =  mfds 
X  (i)  =  x  (i)  +dxi 

dobj  =  fn  f  (X) 

df  (1)  =  (Hot  j-obj) /dxi 

x(i)  =  xd(i) 
next  1 
return 

rem  calculate  the  DG(i,j) 
gosub  3400 
for  i  =  1  to  ndv 

dxi  =  f  dm*x  (i) 

if  dxi  <  mfds  then  dxi  =  mfds 

X  (i)    =    X fi)  +dxi 

for    j    =    1    tc   nave 
k   =    iwrk  (j) 
deon   =    fn_g(x,k) 
dg(j,i)     =    (dcon-wrkl  ( j)  ) /dxi 

next    1 

x(i)    =    xO(i) 
next    1 
return 

rem   ealcilate    the   push-off    factor 
for    i   =    1    to   nave 

thta(i)     =    thtO*  (1. -wrki  (i) /ace)  **2 

if    thta  (i)     >   thtm    then   thta  (i)    =    thtm 
next    i 
return 

rem   normalize    the    DF(i) 
gosub    3200 
rsg   =    0. 
for    i    =    1    to   ndv 

fsq   =    fsg  +  df  (i)  **2 
next    1 

fsq  =  sqr  (fsc) 
if  fsc  =  0.  then  fsg  =  zro 
for  i  =  1  to  ndv 

wrk3(i)  =  (1./fsq)*df  (i) 
next  i 
return 

rem  normalize  the  DG (i) 
gosub  3250 
for  i  =  1  to  nave 

gsq  =  0 . 


46 


U21 5     for  j  =  1  tc  ndv 

4220        gsq  =  gsq+dg (i, j) **2 

U225     next  j 

4230     gsq  =  sqr  (gsq) 

4232     if  gsq  =  0.  tnen  gsq  =  zro 

4235     for  j  =  1  to  ndv 

4240       wrky{i,j)  =  ( 1. /gsq)  *dg  (i,  j) 

4245     next  j 

4250  next  i 

4255  return 

4400  rem  exist  the  violate  constraints 

4405  gosub  3350 

4410  for  i  =  1  to  nave 


4420     for  i  =  1  tc  ndv 
4430       a 
4440     next 


1  =  T  tc  nav 
4430       a"ti,j)  =  •wrky(i,J) 


44 4 u     next  1 

4450     a(i,ndv  +  1)  =  thta(i) 

4460  next  i 

4470  for  i  =  1  to  ndv 

4480     p  (i^  =  -wrk3  (i) 

4490  next  i' 

4500  £(ndv+1)  =  phid 

4510  for  i  =  1  to  nave 

4520    yy  =  0 

4530     for  j  =  1  tc  ndv+1 
4540        XX  =  a(i,j)  *p  ( j) 

4560     next  i 

4570   c(i)  =  (-i.)*yy 

4580  next  i 

4536  ndt  =  nave 

4590  return 

4600  rem  only  exist  aetive  constraints 

4605  gosub  3350 

4610  for  i  =  1  to  nave 

4620    for  i  =  1  tc  ndv 

4630       a1i,j)  =  wrky(i,j) 

4640    next  i 

4650     a(i, ndv  +  1)  =thta(i) 

4660  next  i 

4670  for  j  =  1  to  ndv 

4680    a(navc+1,j)  =  wrk3(j) 

469  0  next  j 

4700  a  (navc+1 ,ndv+1)  =  1. 

4710  p(ndv+1)  =  1. 

4720  for  i  =  1  to  nave+1 

4730    ce  =  a  (i, ndv  +  1) *p (ndv  +  1) 

474  0    e(i)  =  (-1.)  *ce 

4750  next  i 

4760  ndt  =  navc+1 

4770  return 

5000  rem   calculate  the  usable-feasible  direction 

5002  gosub  3000 

5005  gosub  3250 

5010  gosub  3450 

5040  for  i  =  1  to  ndt 

5050     for  j  =  1  to  ndv+1 

§^^9        y^Hy  (D/i)  =  a  (i,  j) 

5070    next  i 

5080  next  i 

5090  for  i  =  1  to  ndt 

5100     for  1  =  1  to  ndb 

5110       ff  =  0. 

5120       for  k  =  1  to  ndv+1 

5130  tf  =  a(i,k)*wrky  (k,1) 

5140  ff  =  ff+{f         -^ 

5150       next  k 

5160       b(i,j)  =  (-1.)*ff 

5170     next  j 
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5180 
5190 
5200 
5210 
5220 
5230 

52ao 

5250 
5260 
5270 
5280 
5290 
5300 
5310 
5320 
5330 
53U0 
5350 
5360 
5370 
5380 
5385 
5400 
5405 
5410 
5420 
54  3  0 
5440 
5460 
5470 
5480 
5490 
5500 
5505 
5510 
5520 
5525 
5530 
5540 
5550 
5560 
5565 
5570 
5580 
5600 
5610 
5620 
5630 
5640 
5650 
5660 
5670 
5680 
5690 
5700 
5710 
5720 
5730 
5740 
5750 
5755 
5760 
5770 
5780 
5820 
6000 
6005 
6010 
6015 


next  i 
iter  =  0 
nmax  =  5*ndb 
rem  begin  iteration 
iter  =  iter+1 
cbmx  =  0 . 
ichk  =  0 

for  i  =  1  to  ndt 
ci  =  c (i) 
bii  =  b  (i,  i) 
if  bii  =  6.  then  5340 
if  ci  >  0.  then  5340 

cb  =  ci/bii 
if  cb  <=  cbmx  then  5340 
ichk  =  i 
cbmx  =  cb 
next  i 

if  cbmx  <  zro  or  iter  >  nmax  then  5550 
if  ichk  =  0  then  5550 

j j  =  iwrk  (ichk) 
if  j1  =  0  then  iwrk  (ichk)  =  ichk  else  iwrk  (ichk)  =  0 
11  b (ichk- ichk)  ,=  0.  then  b  (ichk, ichk)  =  zro 
bb  =  l./b  (ichk, ichk) 
if  bb  =  0.  then  bb  =  zro 
for  i  =  1  to  ndb 

b  (ichk,i)  =  bb*b(ichk,i) 
next  i 

c (ichk)  =  cbmx 
for  i  =  1  to  ndb 

if  i  =  ichk  then  5530 

bbi  =  t(i,ichk) 

b(i,ichk)  =  0. 

for  j  =  1  to  ndt 

if  1  =  ichk  then  5520 

6(i,j)  =  b(i,j)-bbi*b(ichk,j) 
next  1 
c  (i)  =  c  (i) -bbi*cbmx 
next  i 
goto  5220 
ner  =  0 

for  i  =  1  to  ndb 
u  (i)  =  0. 
i  =  iwrk  (i) 
if  1  >  0  then  u  (i)  =  c  (j) 
next  i 

for  i  =  1  to  ndt 
ff  =  0. 
for  1  =  1  to  ndb 

ff  =  ff+wrky  (i,  j)*u(  j) 
next  j 

next  1 

return 

reir  normalized  the  s(i) 

ssq  =  0. 

for  i  =  1  to  ndv 

ssq  =  ssg  +  s  (i)  **2 
next  1 

ssg  =  sgr  (ssg) 

if  fslp  =  0.  then  fslp  =  zro 
for  i  =  1  to  ndv 

s(i)  =  (l./ssg)  *s(i) 
next  i 
return 

rem   one-dimensicnal    search    for    initial   feasible    point 
rem   calculate    for    slope    of    f (x) 
fslp    =   0. 
for   i   =    1    to    ndv 
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6020 
6025 
6035 
604  0 
6045 
6050 
6055 
6060 
6065 
6067 
6070 
6075 
6076 
6080 
6085 
6090 
6095 
6100 
6105 
6110 
6115 
6120 
6125 
6130 
6135 
6140 
6145 
6150 
6155 
6160 
6165 
6170 
617^ 
6180 
6185 
6190 
6195 

6200 
6205 
6210 
6215 
6220 
6225 
6230 
62  3  5 
6237 
6240 
6245 
6250 
6252 
6255 
6260 
6265 
6270 
6275 
6280 
6285 
6290 
6295 
6300 
6305 
6310 
6315 


6  32  1 
6325 
6326 


fslp 
next  i 
rem  iden 
a  lev  =  0 
flow  =  o 
for  i  = 

wrkl  { 
next  i 
rem  find 
if  fslp 
a 1st  =  a 
for  i  = 

if  s( 

walp 

if  wa 
a1 
next  i 
rem  upda 
alph  =  a 
gosub  70 
gosub  71 
rem  calc 
fist  =  f 
for  i  = 

wrkl  ( 
next  i 
rem  chec 
ncvl  =  0 
for  i  = 

if  wr 
nc 
next  i 
if  ncvl 
alst  =  0 
goto  610 
rem  find 
rem  2-po 

for 
aO  =  flo 
a1  =  fsl 
a2  =  (f1 
if  a2  <= 

a2nd 
rem  2-po 
for  i 

aO 

if 

a1 

if  a1 

walp 

walp 

if  wa 
a2 
next  i 
rem  upda 
alph  =  a 
gosub  70 
gosub  71 
rem  calc 
f2nd  =  f 
for  i  = 

wrk2  { 
next  i 
rem  find 
3-po 
f1  =  flo 
alpl  =  a 
f2  =  fis 
alp2  =  a 


a1 


if 


=  fslp+df (i) *s  (i) 

fy  the  initial  point. 

1  to  nice 
i)  =  gcv(i) 

alst    •   the    1st    mid-point. 
=    0.     then    fslp    =    zro 
boi*f  low/abs  (fslp) 
1    to    ndv 

i)    =    0 .    then    s  (i)    =   zro 
=    alpx*x  (i)/abs  fsji)  ) 
Ip   >    alst   then   60^5 
St   =    walp 

te   x   for   alst. 

1st 

00 

00 

ulate    the    fist    and    wrkl(i) 

n    f(x)  . 

1    to   nigc 

i)    =    fn_g(x,i) 

k    the  -feasibility. 

1    to    niqc 

k1  (i)     <   vcc    then   6  170 

v1    =    ncv1+1 

=    0    then    6200 

.5*a1st 

5 

a2nd    ;    the    2nl    mid-point, 
ints    quadratic    fit    interpolation 
minimum  f  (alpa)  . 
w 

St-a1*a1st-a0) /(a1st**2) 

0.    then   a2   =    zro 
=    -a1/  (2.*a2) 

ints    linear   interpolation    for    g(alpa)=0 
1    to   nice 
wrkl  (i) 

St   =    0.    then    alst    =    zro 
(wrkl  (i)-aO)/a1st 

<=    0.    then   a1    =   zro 
=    -aO/al 

<=    0.     then    walp    =    1000. 
Ip    >=    a2nd   then    6265 
na   =    walp 

te   X    for   a2nd. 

2nd 

00 

00 

ulate    f2nd   and    wrk2  (i) 

n  f(x). 

1  to  nigc 

i)  =  fn_g(x,i) 

final  point  aupr  by  using 
ints  quadratic  fit. 
w 

low 
t 
1st 
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6330 
6331 
63  3  5 

63ao 

6342 
6345 
6347 
6350 
6355 
6360 
6  36  5 
6370 
637  5 
6376 
6377 
6380 
6385 
6390 
6395 
6400 
6405 
6410 
6415 
6420 
6425 
6430 
6435 
6440 
6445 
6450 
6455 
6460 
6465 
6470 
6475 
6480 
6485 
6490 
6495 
6500 
6505 
6510 
6515 
6520 
6525 
6530 
6535 
6540 
6545 
6550 
6555 
6560 
6565 
6567 
6569 
6571 
6573 
6575 
6577 
6579 
6600 
6603 


f3  =  f2n 
alp3  =  a 
rem  3-po 
gosub  66 
if  a2  = 
a3rd  =  - 
if  a3rd 
for  i  = 
f1  = 
f2  = 
f3  = 
gosub 
gosub 
if  alps 

a3 
next  i 
rem  upda 
alph  =  a 
gosub  70 
gosub  71 
rem  calc 
fupr  =  f 
for  i  = 

wrku  ( 
next  i 
rem  find 
f1  =  fis 
f2  =  f2n 
f3  =  f3r 
alpl  =  a 


d 
2nd 


ints    quadratic    fit    interpolation 


1000. 


alp2  = 
alp3  = 
rem  3- 
gosub 
if 


a 

a 

po 

56 

a2 


aupr 
for  1  = 
f1  = 
f2  = 
f3  = 
alpl 
alp2  =  a 
alp3 
gosub 
gosub 
if  al 
aupr 
next  1 
rem  upda 
alph  =  a 
gosub  70 
gosub  71 
rem  eval 
fupr  =  f 
for  i  = 

vrku 
next  i 
reir  find 
gosub  14 
return 
rem  3-po 
if  alpl 
then  ret 


6605  a2  = 


6610 
6615 
6620 
6630 


return 
rem  zero 


00 

0.  then  a2  =  zro 

al/f  2.*a2) 

<=  0.  then  a3rd  = 

1  to  niqc 

wrklfi* 

wrki (i 

wrk2  (i 

6600 

6630 
>  a3rd  then  6380 
rd  =  alp^ 

te  x  for  aupr 

3rd 

00 

00 

ulate  the  fupr  and  wrku(i) 

n  f(x)  . 

1  to  nigc 

i)  =  fn_g  (x,i) 

4th  new  point, 
t 
d 
d 

1st 
2nd 
3rd 

ints  quadratic  fit. 
00 

=  0.  then  a2  =  zro 
=  -a1/(2.*a2) 
1  to  niqc 
wrki  (i) 
wr k2  (i 
wrk3  (i' 
=  alst 
2nd 
=  a3rd 

6600 

6630 
ps    >    aupr   then    6540 
=    alps 

te  x    for   aupr 
upr 

o6 

00 

uate    fupr   and    wrku(i> 

nf(x). 

1    to    niqc 

(i)    =    fn_g(x,i) 

optimum  alpa  for  not  violating  constraints 
306 

ints  quadratic  fit. 

=  alp2  cr  alp2  =  alp3  or  alpl  =  alp3 

urn 

3-f  1)/  (alp3-alp1)- 

-f  1)  /(alp2-alp1)  )/ (alp3-alp2) 

-f  1)/  alp2-alp1)  -a2*(alp1+alp2) 

a1*alp1-a2*alp1**2 

of  polynomial  for  g(alpa) 
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6635  dd  =  a1**2-U.*a2*aO 

6640  if  dd  <  0.  then  6715 

6642  if  a2  <=  0.  then  a2  =  zro 

6645  if  a2  =  0.  then  a2  =  zro  ^ 

6650    alpb  =  (-a1 +£gr  (dd)  )  /  (2.  «a2) 

6655     alpc  =  (-al-sgr  idd)  j  /  (2.  *a2) 

6660     if  alpb  <=  0  and  alpc  <=  0.  then  6715 

6665     if  alpb  >=  0.  and  alpc  >=  0.  then  6695 

6670     if  alpb  >=  0.  and  alpc  <  0.  then  6685 

6675       alps  =  alpc 

6680  goto  672D 

6685       alps  =  alpb 

6690     goto  5720 

6695     if  alpb  >=  alpc  then  6710 

6700       alps  =  alpb 

6705     goto  6720 

6710       alps  =  alpc 

6712     goto  6720 

6715  alps  =  1000. 

6720  return 

6780  rem  update  aboj  and  alpx 

6790  delf  =  abs  (obj-nobi) 

6795  diff  =  abs  (delf/obji 

6800  abcj  =  fabo j+dif f ) /2. 

6815  walp  =  0. 

6816  welx  =  0. 

6820  for  i  =  1  to  ndv 

6830     delx  =  abs  (xO (i) -x (i)  ) 

6850     difx  =  abs(delx/xO (i)  ) 

6855    if  delx  >=  welx  then  welx  =  delx 

6860     if  difx  <=  walp  then  6880 

6870       walp  =  difx 

6880  next  i 

6890  alpx  =  (alpx+walp) /2. 

6910  dabf  =  accf *abs  (ob j) 

69  90    return 

7000    rem   update   the    x(i) 

7010    for  i   =    1    to   ndv 

7020  X  (i)    =    xO  (i)+alph*s  (i) 

7030    next    i 

7040    return 

7100    rem   check    the    side-constraints. 

71 10    for   i   =    1    to   ndv 

7120     if  x(i)  <  lcwb(i)  then  x(i)  =  lowb  (i) 

7130     if  x(i)  >  uprb  (i)  then  x  (i)  =  uprb  (i) 

7140  next  i 

7150  return 

8000  rem  estimate  the  alpa 

8010  fstr  =  flow 

8020  alpa  =  alow 

8030  nvcl  =  0 

8040  for  i  =  1  to  nice 

8050     if  wrki  (i)  <  vcc  then  8070 

8060     nvcl  =  nvc1+1 

8070  next  i 

8080  if  nvcl  >  0  then  8120 

8090  if  fist  >  fstr  then  8120 

8100     alpa  =  alst 

8110     fstr  =  fist 

8120  nvcl  =  0 

8130  for  i  =  1  to  niqc 

8140     if  wrk2(i)  <  vcc  then  8160 

8150     nvcl  =  nvc1+1 

8160  next  i 

8170  if  nvcl  >  0  then  8210 

8180  if  f2nd  >  fstr  then  8210 

8190     alpa  =  a2na 

8200     fstr  =  f2nd 

8210  nvcl  =  0 
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8220  for   i  =  1  to  rice 

R230     if  wrk3  (i)  <  vcc  then  3250 

8240     nvcl  =  nvc1+1 

8250  next  i 

8260  if  nvcl  >  0  then  8300 

8270  if  f3rd  >  fstr  then  8300 

8280     alpa  =  a3ra 

8290     fstr  =  f3rd 

8300  nvcl  =  0 

8310  for  i  =  1  to  niac 

8320    if  wrku(i)  <  vcc  then  8340 

8330     nvcl  =  nvc1+1 

8340  next  i 

8350  if  nvcl  >  0  then  8390 

8360  if  fupr  >  fstr  then  8390 

8370     alpa  =  aupr 

8380     fstr  =  fupr 

8390  alph  =  alpa 

8400  return 

9000  rem  one-dimensional  search  for  initial 

infeasible  point, 
9002  ii  =  1 
9004  gcvm  =  wrkl  (1) 
9006  ror  i  =  1  to  nave 
9008  if  wrk1(i)  <=  gcvm  then  9014 
9010     ii  =  i 
9012     gcvm  =  wrkl  (i) 
9014  next  i 

9016  rem  calculate  the  slope  of  badly  violation. 
9018  gslp  =0. 
9020  for  i  =  1  to  ndv 
9022     gslp  =  gslp+dg(ii,i) *s(i) 
9024  next  i 

9026  rem  calculate  the  alph. 

9027  if  gslp  =  0.  then  gslp  =  zro 

9028  alph  =  -gcvm/gslp 
9030  rem  update  X  for  alph. 
9032  gosub  7000 

9034  gosub  7100 

9036  rem  evalute  the  objective  and  constraint. 

9038  obi  =  fn  f  (x) 

9040  for  i  =  T  to  niqc 

9042    gcv  (i)  =  fn  g  (x,i) 

9044  next  i        -^ ^  '  / 

9046  rem  calculate  the  NVC. 

9048  gosub  3500 

9050  if  nvc  =  0  then  return 

9052  rem  update  initial  value. 

9054  for  i  =  1  to  ndv 

9056     xO  (i)  =  x  (i) 

9058  next  i 

9060  rem  calculate  df (i) ,dg  (i, j)  and  push-off  factor. 

9062  gosub  3800 

9064  gosub  3900 

9066  gosub  4000 

9068  rem  normalize  the  df (i) , dg  (i, 1) 

9070  gosub  4100 

9072  gosub  4200 

9074  rem  find  the  search  direction. 

9076  gosub  4400 

9078  gosub  5000 

9030  goto  9000 

9200  rem  print  the  results 

9205  print  '» 

9210  print  '********«**  data  ************ 

9215  print  " 

9220  print  'The  number  of  design  variables      =  '.ndv 

9225  print  'The  number  of  inequality  constraints  =  ',niqc 

9230  print  " 
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9235 
9240 
9245 
9250 
9255 
9260 
9265 
9270 
9275 
9280 
9285 
9290 
9295 
9  30  0 
9305 
9310 
9315 
9500 
9510 
9520 
9530 
9540 
9550 
9560 
9570 
9580 

9590 
9600 
9610 
9620 
9630 
9640 
9650 

9660 
9670 
9680 
9690 

9700 
9800 


print 
print 
print 
for  i 

next  1 

print 

print 

print 

print 

print 

print 

print 

tor  i 

pri 
next  i 
return 
rem  de 
mxit  = 
fdm  = 
mf  ds  = 
vcc  = 
ace  = 
thtO  = 
thtm  = 
phid  = 

accf  = 
accx  = 
zro  = 
espl  = 
bnlo  = 
bnup  = 
dalp  = 


1 1 


]The  objective  value  =  ',obj 

ables  ****** 

(i) 


******  design  vari, 

=  1  to  ndv 

nt  •  x('  ;i;  ')  =  '  ,x 


'the  number  of  active  constraints  = 

'the  number  of  vio 

*****  constraint  v 

=    1    to    niac 

nt    *g('  ;i; ')    =    '  ;gcv  (i) 


;  nac 


Date   constraints   =    '  ;nvc 
alue   ****' 


1 


fault 

50 
.01 

.00 
.004 
-.  1 

1. 

50. 

10000 

.001 
0.001 
0001 
.005 
1.e+ 
e+7 
1 


abcj  = 
alpx  = 
ndvm  = 
nigra  = 


return 
end 


=  1 

=  .0 


=  0.  1 


1 
21 
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nunber 
maxiro 
finit 
maxim 
viola 
activ 
push- 
maxim 
!  we 
wh 
1  abs 
!  abs 
!  def 
used 
!  th 
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